#include "MultigridSolver.h"
#include "FEMSpace.h"
#include <math.h>
#include <ctime>
#include <fstream>
#include <iostream>
#include <string>
#define pi 4.0*atan(1.0)

// double u(double *p)
// {
//     return sin(pi*p[0])*sin(pi*p[1]);
// }
// double f(double *p)
// {
//     return 2 * pi*  pi* sin(pi * p[0]) * sin(pi * p[1]);
// }

// double bc(double *p)
// {
//     return sin(pi*p[0])*sin(pi*p[1]);
// }

double u(double *p)
{
    return sin(p[0])*sin(p[1]);
}
double f(double *p)
{
    return 2*sin(p[0]) * sin(p[1]);
}

double bc(double *p)
{
    return sin(p[0])*sin(p[1]);
}


int main(int argc, char* argv[])
{
    //g++ -o main test_solver.cpp -std=c++11 -I /usr/include/eigen3/ -g
    int n,smoothetype,order;
    std::string elementtype;
    std::cin>> n >> smoothetype >> elementtype >> order;
    RectangleDomain* r = new RectangleDomain({{0,0},{1,0},{1,1},{0,1}});
    std::vector<Boundary<2> > B = r->boundary();
    Mesh<2> *m;
    Element<2> *e;
    if(elementtype == "P"&& order ==1)
    {
	m = new P1Mesh(r, {POW2(n), POW2(n)});
	e = new Triangular_1_Element();
    }
    else if(elementtype == "P"&& order ==2)
    {
	m = new P2Mesh(r, {POW2(n), POW2(n)});
	e = new Triangular_2_Element();
    }
    else if(elementtype == "Q"&& order ==1)
    {
	m = new Q1Mesh(r, {POW2(n), POW2(n)});
	e = new Quadrilateral_1_Element();
    }
    else if(elementtype == "Q"&& order ==2)
    {
	m = new Q2Mesh(r, {POW2(n), POW2(n)});
	e = new Quadrilateral_2_Element();
    }
    
    //Mesh<2>* m = new P1Mesh(r,{POW2(n), POW2(n)});
    //Element<2>* e = new Triangular_1_Element();
    //Mesh<2>* m = new Q1Mesh(r,{POW2(n), POW2(n)});
    //Element<2>* e = new Quadrilateral_1_Element();
    //Mesh<2>* m = new Q2Mesh(r,{POW2(n), POW2(n)});
    //Element<2>* e = new Quadrilateral_2_Element();
    //Mesh<2>* m = new P2Mesh(r,{POW2(n), POW2(n)});
    //Element<2>* e = new Triangular_2_Element();
     PoissonEquation* equ = new PoissonEquation();
    equ->SetRightHandsTermFunction(f);
    equ->SetBoundaryConditionFunction(bc);
    BoundaryFunction<2>* dir = new Dirichlet<2>(bc,B);
    BoundaryCondition<2> bfc(dir);
    Poisson_2D poissonproblem(m, e, equ, bfc);
    poissonproblem.AssembleStiffMatrix();
    poissonproblem.AssembleRightHandsTerm();
    //possionproblem.DealWithBoubdaryCondition();

    MultigridSolver* Solver;
    if(elementtype == "P"&& order ==1)
    {
	Solver =new P1_MGSolver(poissonproblem.A(), poissonproblem.Rhs(), bc, m, n);
    }
    else if(elementtype == "P"&& order ==2)
    {
	Solver = new P2_MGSolver(poissonproblem.A(), poissonproblem.Rhs(), bc, m, n);
    }
    else if(elementtype == "Q"&& order ==1)
    {
	    Solver =new Q1_MGSolver(poissonproblem.A(), poissonproblem.Rhs(), bc, m, n);
    }
    else if(elementtype == "Q"&& order ==2)
    {
	Solver = new  Q2_MGSolver(poissonproblem.A(), poissonproblem.Rhs(), bc, m, n);
    }
    
    Solver->setTolerance(1e-16);
    Solver->set_smoothetype(smoothetype);
    Solver->compute(poissonproblem.A());
    Eigen::VectorXd real_solution = Eigen::VectorXd::Zero(m->n_dofs());
    std::cout << m->n_dofs() << std::endl;
    double point[2];
    for(int j = 0;j < m->n_dofs();j++)
    {
	Dofs<2> dof = m->DofsofIndex(j);
	point[0] = dof[0];
	point[1] = dof[1];
	real_solution(j) =  u(point);
    }
    //std::cout << real_solution <<std::endl;
     Eigen::VectorXd solution = Solver->Solve(real_solution);
     std::ofstream output("./output/result2.txt");
     for(int j=0; j <solution.size();j++)
      {
    	 Dofs<2> dof = m->DofsofIndex(j);
    	 output << dof[0] << " "<<dof[1] << " " << solution[j] << std::endl;
     }
     output.close();
    return 0;
}
    
